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Abstract 

This paper presents numerical simulations for a 
high-speed transonic compressor rotor, NASA Rotor 
35, using an unstructured grid flow solver U 2 NCLE. 

It solves Reynolds averaged Naiver-Stokes equations 
using an arbitrary Mach number solution algorithm. 
Experiments were conducted at the NASA Gleen Re- 
search Center to investigate compressor stability en- 
hancement using steady and discrete tip injections. 
The mass injection upstream of the tip of a high- 
speed axial compressor rotor has shown effectiveness 
in suppressing stall in tip-critical rotors. This work 
presents numerical simulations of a high-speed com- 
pressor rotor with and without tip injections. Methods 
to model the flow injection are described. Numerical 
simulations are conducted at 80 and 100 percent of 
the design speed. Results show that numerical solu- 
tions are consistent with the experimental trends. 

Introduction 

A research effort supported by the NASA Glenn 
Research Center was initiated in 2001 for the purpose of 
exploring the use of an unstructured solver in die cal- 
culation of turbomachinery flows, and specifically, to 
consider the possible benefits of developing an unstruc- 
tured version of the MSU TURBO code [1]. Whereas 
the present MSU TURBO code, distributed to U.S. en- 
gine companies by NASA Glenn, is a heavily used par- 
allel unsteady compressible Reynolds averaged Navi- 
er-Stokes flow solver, it is a multiblock code that uses 
structured subgrids. Consequently, there are some geo- 
metric configurations that are difficult and time con- 
suming to grid. For example, the Army Research Labo- 
ratories Vehicle Technology Directorate (VTD) has an 
urgent need for the capability to numerically simulate 
the flow through centrifugal compressors with vaned 
diffusers, which are difficult to grid using structured 
grid technology due to the complexity of the geometry. 
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In the past few years an unstructured Reynolds averaged 
Navier-Stokes flow solver, referred to as U 2 NCLE (Un- 
structured UNsteady Computation of fieLd Equations), 
has been developed in the Computational Simulation 
and Design Center at Mississippi State University for 
compressible, incompressible, and arbitrary Mach num- 
ber flow regions [2-4]. It is prudent to consider the use 
of this new technology for application to turbomachin- 
ery. 

The test case studied in this work is a high-speed 
axial compressor rotor, NASA Rotor 35. Experiments 
to enhance the stability of an axial compressor rotor us- 
ing tip injection techniques were conducted by Suder, et 
al. [5] at the NASA Glenn Research Center. They found 
that mass injection upstream of the tip of a compressor 
rotor was effective in suppressing stall in tip-critical ro- 
tors. Experiments also showed that stability enhance- 
ment was related to the mass-averaged axial velocity at 
the tip, and the configuration of the injector. To investi- 
gate the impact of injector parameters on stability im- 
provement for an axial compressor rotor, Suder, et al. [5] 
performed numerical simulations using the Average 
Passage code (APNASA) developed by Adamczyk [6]. 
Time-accurate unsteady simulations for a transonic 
compressor stage. Stage 35, with tip injections were re- 
cently conducted at the Computational Simulation and 
Design Center at Mississippi State University using the 
MSU-TURBO code [7]. These investigations indicate 
that the effect of the injection is sensitive to the momen- 
tum and velocity of the flow leaving the injector. The 
injector effectiveness is maximized when the injected 
flow is choked. 

The purpose of this work is to assess the capability 
and accuracy of the current unstructured grid flow solv- 
er U 2 NCLE to compute flows in turbomachinery. It is 
intended to develop a methodology to simulate flow in- 
jections and analyze the effect of stability enhancement 
on high-speed compressor rotors and stages. The 
technology developed in this work will later be used to 
investigate the stall inception process and stall control 
techniques in a high-speed centrifugal compressor 
stage being tested at the US Army Research Laboratory. 
To this end, two computational meshes are generated for 
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simulations of an isolated Rotor 35 with and without tip 
injections. There are 12 injectors uniformly spaced 
around the compressor annulus upstream of the rotor tip 
in the experiment. However, this work solves only one 
injector and three blade passages based on the periodici- 
ty of the flow field. Numerical simulations are con- 
ducted at 80 and 100 percent of the design speed. All 
computational results, including with and without injec- 
tion flows, are verified with the experimental data. 
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Numerical Methods 


To simulate complex viscous flows in turboma- 
chinery, a parallel version of an unstructured Reynolds 
averaged Navier-Stokes flow solver U 2 NCLE was de- 
veloped at the Computational Simulation and Design 
Center at Mississippi State University [4]. The solver 
solves preconditioned three-dimensional compressible 
governing equations in absolute reference frame and/or 
rotating reference frame. The numerical scheme, re- 
ferred to as an arbitrary Mach number solution algo- 
rithm [8], allows calculations of both low-speed and 
high-speed flows in turbomachinery. Validation and 
verification of the U 2 NCLE solver for isolated rotors for 
low Mach number and transonic flows have been re- 
ported in [4] [9]. 

Primitive-Variables Conservation Laws 


where E c = (y — i)M 2 r is an Eckardt number. 0 is the 
normal velocity with respect to the surface of a control 
volume. The viscous shear stresses are given as 
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The unsteady three-dimensional compressible Re- 
ynolds-averaged Navier-Stokes equations are cast in a 
rotating Cartesian coordinate system with a rotational 
speed of Q. The conservative flux formulation is writ- 
ten in terms of primitive variables in the absolute frame. 
In order to solve arbitrary Mach number flows, a pre- 
conditioning matrix r ~ 1 is introduced into the primitive 
variable conservation laws, which can be expressed as 


S is the source term representing both body forces S b f in 
rotating coordinates, as well as injected flows S in j y 
which will be described next. 

Note that the above equations are non-dimen- 
sionalized by the following reference values: density, 
g r \ velocity, U r \ temperature, T r , pressure, length, 
U\ time, Lr/U r \ energy and enthalpy, C p T r . 


Mr~' JU qdV+i F . ndA = SdV (D 
J Q J dQ * Q 

where r ~ 1 is a constant diagonal matrix that only de- 
pends on the reference Mach number M r [8] [10] 

r-' = diag [ 1 , 1 , 1 , 1 , l/flM,)] ( 2 ) 

and M is a transformation matrix from conservative 
variables Q = (p, pu , pv, pw, pe ( ) T to primitive vari- 
ables q = (p, u, v, w, p) r , where u, v, and w, are the 
absolute velocity components in jc, y, and z directions, 
respectively, p is the density; p is the pressure; and e t 
is the total energy, n = (n x% n yt n z ) T is the outward 
pointing unit normal to the surface of the control volume 


Source Term 

The source term in the above Eq. (1) has two parts 
S = S h j + 
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where S b f is the body force (Coriolis and centrifugal 
forces) in rotating coordinates, and Q x , Qy> and Q z are the 
components of the rotational speed Q of the coordinate 
system. The second term S in j is used to model the effect 
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of injected flows by adding the extra mass, momentum, 
and energy into the governing equations. m inJ is the in- 
jected massflow rate, u m] , v my , w inj are the injected veloc- 
ity components in x, y, z directions, and e inj is the total 
energy of injected flow. V is the total volume of cells 
that injected flow is added on. An uniform velocity dis- 
tribution is assumed on the injector area. Noting that Sy 
is added to all control volumes in rotating coordinates, 
while S in j is only added to regions that are marked as the 
injector area. 

Numerical Schemes 

Preconditioned equations (1) can be expressed in a 
differential form as 

A + Mr q -'\r q M~'AM\ ( £ = 5 (5) 

where F q M~ x AM is the system matrix for precondi- 
tioned equations, and A = dF/dQ is the flux Jacobian 
with respect to Q. The inviscid flux approximation at 
the face of a control volume for the preconditioned sys- 
tem is 

Ft = \(Hq L ) + F{q R )) 

-±Mr;'\r 9 M~'AM\(q R -q L ) ( 6 ) 

The eigensystem of the above Eq. (6) is evaluated based 
on the averaged variables q. 
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Since the eigensystem and the flux are evaluated based 
on the above averaged variables instead of Roe vari- 
ables, the current approach is considered to be an exten- 
sion of the Roe flux approximation [10]. For first-order 
accurate differencing, quantities q L and q R are set equal 
to the data at the nodes lying on either side of the face. 
For a second-order scheme, these values are computed 
with a Taylor series expansion about the central node of 
the control volume with an unweighted least-squares 
procedure [11]. 

There are two issues in reconstructing the viscous 
scheme. One is the positivity of the discrete operator, 
which requires that all of the stencil weights be positive. 
Positivity is a key property for numerical stability. 
Another key property is the linearity-preserving of gra- 


dients, which means that for a linear distribution of ve- 
locity, one should obtain a constant gradients of veloc- 
ity. In this work, a positive scheme proposed in [8] is 
used to compute gradients of flow variables at the me- 
dial dual faces for stable and accurate simulations of vis- 
cous flows. 

The numerical algorithm is based on a node-cen- 
tered, finite volume implicit scheme on unstructured 
mixed element grids. The time-advancement algorithm 
is based on Newton’s method, which yields a linear sys- 
tem of equations for the solution at each time step. The 
solution of the sparse system of equations is obtained by 
symmetric Gauss-Seidel relaxations. Normally 8~15 
symmetric Gauss-Seidel subiterations are adequate at 
each time step. 

Turbulence Model 

The current U 2 NCLE solver has several turbulence 
model options, including Spalart and Allmaras one- 
equation model [12], and q-Q), k-e, k-a) two-equation 
models [13]. These turbulence models are loosely 
coupled with the mean flow equations, meaning that 
they are solved separately following the mean flow 
solution at each time step. This has an advantage of eas- 
ily modifying the existing models or adding a new tur- 
bulence model to the solver. For the current study, one- 
equation turbulence model of Spalart and Allmaras is 
used in all calculations. Since there is no wall-function 
model used in the solver, grid points need to be packed 
near viscous surfaces in order to fully resolve the bound- 
ary layer. An ongoing effort is to explore the use of two- 
equation models for future turbomachinery simulations. 

It should be noticed that a source term for the eddy 
viscosity is added to the Spalart and Allmaras turbu- 
lence equation to model the flow injection. This is nec- 
essary to maintain a stable solution for the strong mixing 
between the injected flow and the primary flow. The 
amount being added to the turbulence equation is esti- 
mated based on the injected massflow rate in the current 
work, although a more accurate value may be obtained 
from the measured data. 

Boundary Conditions 

The boundary conditions required for turbomachin- 
ery flow computations were derived in detail in Ref [9] . 
These boundary conditions are based on characteristic 
relations that propagate in the direction of characteristic 
waves, i.e. the eigenvalues of preconditioned system of 
equations. Under the current grid topology, all bound- 
ary faces have positive outward pointing normals. A 
positive eigenvalue means that the corresponding char- 
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acteristic wave propagates from the interior domain to 
the boundary, while a negative eigenvalue means that 
the characteristic wave propagates from the outside of 
the computational domain into the boundary. Those 
characteristic variables with negative eigenvalues are 
replaced by the total pressure and total temperature at 
the inflow, and the static pressure or radial equilibrium 
equation at the outflow for turbomachinery flow com- 
putations. The current unstructured flow solver has the 
following boundary conditions related to turbomachin- 
ery simulations: 

• Inlet specifying total pressure, total temperature, 
and flow direction. 

• Exit specifying constant back pressure. 

• Exit with radial equilibrium condition. 

• No-slip adiabatic wall. 

• No-slip constant temperature wall. 

• Rotating viscous surface for shroud or hub. 

• Periodical condition for axisymmetry surfaces 
(require axisymmetry surface grids). 

• Sliding interface for rotor/stator interaction. 

• Source term for modeling injected flows. 

Results 

The test case studied in this work is a high-speed 
transonic compressor rotor, NASA Rotor 35. It has 36 
blades, an inlet tip radius of 25.4 cm, a hub-tip radius 
ratio of 0.7, and aspect ratio of 1 . 19, a tip solidity of 1 .3, 
and an axial chord of 2.72 cm at the tip and 4.12 cm at 
the hub. The compressor was designed for axial inlet 
flow, and inlet relative Mach number is 1.48 at the tip 
at the design speed. The rotor tip clearance at design 
speed is 0.86 percent of the tip chord. The compressor 
aerodynamic design and blade coordinates were re- 
ported by Reid and Moore [14]. 

Computations were performed at 80 and 100 per- 
cent of the design speed. The reference values at the in- 
let of the rotor are used to normalize the flow variables 
in the code. The reference Mach number is 0.52. The 
reference length is the rotor diameter. The Reynolds 
number is 5.2434x10* based on the entrance velocity 
and the rotor diameter. Standard day pressure and tem- 
perature are used as the total conditions to specify the 
inflow boundary condition. Various back pressures are 
applied at the exit to obtain different operating condi- 
tions. All computations were carried out in parallel 
mode using Message Passing Interface (MPI) for inter- 
process communication. Simulations were performed 
on a large-scale Linux cluster (1024 Intel Pentium III 
1.3 GHz processors and 5 12 MB RAM per processor) lo- 


cated at the Engineering Research Center at Mississippi 
State University. 

Rotor 35 without Tip Injection 

The first case studied is an isolated rotor without tip 
injection. Simulations were conducted to obtain the 
baseline operating characteristics at 80 percent of the 
design speed (a peak efficiency operating condition). 
This corresponds to a speed of 13,755 rev/min and a ro- 
tor tip speed of 363 m/s. Laser Doppler Velocimeter 
(LDV) measurements [15] were acquired in the NASA 
Glenn single-stage axial flow compressor facility using 
the NASA Rotor 35 operating in a rotor-only configura- 
tion, and are used to verify the computed results. The 
computational grid (Figure 1) was generated for a single 
blade passage using a grid generation tool GUM-B [16]. 
A periodical boundary condition is applied at the axi- 
symmetric surfaces. In order to simulate the tip clear- 
ance flow, a true gap was generated between the rotor tip 
and the casing, see Figure 2. The grid has approximate- 
ly one million points, which is partitioned into 32 block 
for parallel computations. Simulations were conducted 
in the rotating reference frame by solving absolute ve- 
locity components. Solutions were considered con- 
verged when the difference between the inlet and the 
exit massflows was below one percent of the total mass- 
flow. 



Figure 1 Single Passage Grid for Rotor 35 



Figure 2 Cutting Plane Grid in Tip Clearance 
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Rotor 36 Operating Rang* at tt% Spood 



Figure 3 Computed and Measured Rotor 35 Operating 
Range at 80% of Design Speed 


Figure 3 shows the comparison of measured and. 
computed total pressure ratios for Rotor 35 at 80 percent 
of the design speed. The computed operating curve was 
obtained by gradually increasing the back pressure at 
the exit in the computational domain, starting from a 
value of 96,8 1 8 N/m 2 (Pa) on the shroud. The measured 
total pressure ratio and efficiency at peak operating con- 
dition are 1.44 and 92 percent respectively, which corre- 
sponds to a massflow of 17.35 kg/s. The computed total 
pressure ratio and efficiency which are closest to the ex- 
perimental data are 1.46 and 90.4 percent, respectively. 
The computed and measured stall points are marked 
with arrows in Figure 3. It seems that computations 
overpredicted the total pressure ratio at the near stall 
point by 3 percent. The predicted stable operating range 
is smaller than the measured value (the computed stall 
massflow is larger than the experimental data), which 
could be caused by the grid resolution and the tip treat- 
ment. Van Zante, et al. [15] studied three different tip 
grid topologies and their impact on the predicted stabil- 
ity. They found that the trajectory of the clearance vor- 
tex can be influenced by different tip grid treatments, 
which in turn affect the stability of the rotor. However, 
the study on grid resolutions was not conducted in the 
current work, and will be conducted in the future. Nev- 
ertheless, the above computed operating characteristic 
curve matches with the experimental trend. 

The tip clearance vortex has a great impact on the 
rotor performance. Hoying, et al. [17] and Adamczyk, 
et al. [18] suggest that rotor stall occurs when the tip 
clearance vortex spills forward of the leading edge. It 
is very important that the tip vortex and its trajectory are 



Experiment Computation 

Figure 4 Computed and Measured Axial Velocity (m/s) 
on the 92% Span Stream Surface at 
Near Peak Efficiency 



Figure 5 Velocity Vectors on x-r Cutting Plane 

correctly predicted in the numerical simulations. 
Figure 4 shows computed contours for axial velocity for 
the 92 percent span stream surface, and measured con- 
tours by Van Zante, et al. [15]. These plots clearly illus- 
trate the lower extent of the flow field region impacted 
by the tip clearance flow. The axial velocity is chosen 
here since the footprint of the clearance flow shows 
most clearly as gradients in axial velocity on this stream 
surface. Figure 5 shows computed velocity vector on 
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Figure 6 Computed and Measured Exit Variables at 
Peak Efficiency 


may be due to the lack of grid resolution in the region. 
Van Zante, et al. [15] suggested that the induced vortex 
has a strong impact on the rotor stability. The rotor sta- 
bility can be enhanced by increasing the strength of the 
induced vortex. It is believed that a refined grid in the 
tip region may help capture the induced vortex, and thus 
improve the predicted operating characteristics as well. 

Figures 6 and 7 show comparison of computed and 
measured total pressure P„ total temperature T t , static 
pressure p s , adiabatic efficiency rj , absolute tangential 
velocity V„ and axial velocity V x at the exit of the blade 
passage along the spanwise direction. These quantities 
are obtained by a massflow averaging procedure along 
the circumferential direction at each span location. 
Comparisons are made at the peak efficiency and the 
near stall point. Computed total pressures (and static 
pressures) at both operating points overpredict the mea- 
sured values, which is consistent with predicted results 
shown in Figure 3. The overpredicted pressure leads to 
slight deviation of adiabatic efficiency. However, com- 
puted results in general match well with the experimen- 
tal data. 
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Figure 7 Computed and Measured Exit Variables at 
Near Stall Point 

the x-r cutting plane at the middle chord location, with 
the axial flow coming from the right. A tip vortex is 
clearly observed in the numerical solution. The trajec- 
tory of the clearance vortex is determined by the wall- 
bounded shear layer near the shroud. The computed 
clearance vortex matches with the numerical results ob- 
tained by Van Zante, et al. [15]. However, an induced 
vortex captured in Ref. [ 15] is not observed here, which 


Rot or 35 with Tip In j ection 

As mentioned in the introduction, one goal of this 
work is to develop a methodology to simulate flow in- 
jections in high-speed transonic compressors, and to in- 
vestigate the stabilizing effect of tip injections on com- 
pressor performance. Casing-mounted injectors are 
located 55 cm (about 200 percent of rotor tip axial 
chord) upstream of the rotor. There are 12 injectors uni- 
formly spaced around the compressor annulus, which 
are designed to generate a jet along the casing upstream 
of the rotor. Both numerical simulations and experi- 
mental measurements were conducted by Suder, et al. 
[5] to investigate the compressor stability enhancement 
using the tip injection. The impact of tip injection on the 
stability limit of an axial compressor stage was investi- 
gated by Weigl, et al. [19]. Their studies showed that jet 
flows upstream of the rotor tip were effective in sup- 
pressing the stall inception in the axial compressor rotor 
and stage. The measured extension of operating range 
for an isolated rotor was about 30 percent at 70 percent 
of the design speed. The extension of operating range 
at 100 percent of the design speed was less than what 
was achieved at 70 percent of the design speed, but was 
still noticeable. Figure 8 shows impact of the tip Injec- 
tion on the stability limit for a transonic axial compres- 
sor stage as measured by Weigl et al. [19]. 

The current work simulates a steady tip injection 
upstream of an isolated Rotor 35 at the design speed, us- 
ing an unstructured grid flow solver. The standard day 
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Figure 8 Impact of Tip Injection on the Stability 
Limit of a Transonic Axial Compressor 
Stage Measured by Weigh et al. [19] 


corrected speed for the design flow condition is 
17,188.7 rev/min with a rotor tip speed of 454.14 m/s. 
The effective injector size in the grid is 4.73 mm in 
depth and 40.23 mm in width, which covers about 75 
percent of the actual injector exit area. An uniform ve- 
locity across the injector exit is assumed. The injected 
flow is aligned with the inlet annulus flow in the down- 
stream axial direction to minimize mixing between the 
annulus flow and the injected flow. The stagnation pres- 
sure for the injected flow is 1.88 times the standard day 
pressure, and the stagnation temperature of the injected 
flow is the same as the standard day temperature. A 
choke flow is assumed at the injector exit, which results 
in a mass flow of 1.026 kg/s (about 5 percent of rotor 
choke flow). Figure 9 shows the computational mesh 
which includes one injector and three rotor-blade pas- 
sages in a 30-degree segment. There is an interface be- 
tween the injector and the rotor which divides the com- 
putational mesh into two components, the injector grid 
and the rotor grid. The grid points are smoothly con- 
nected between the injector and the rotor. Simulations 
are implemented in a steady-state mode, in which the 
flow in the injector grid is solved in the fixed absolute 
frame, while the flow in the rotor grid is solved in the 
rotating frame. Since the absolute velocity components 
are cast in the governing equations in both absolute and 
rotating frames, fluxes at the interface of the injector 
and the rotor are directly computed. Periodical bound- 
ary conditions are assumed at axisymmetric surfaces for 
both injector and rotor grids. Because there is no rela- 
tive motion between the injector and the rotor, the above 
implementation can be viewed as to solve an isolated 
axial rotor with a stationary inflow source (injector) up- 
stream of the rotor in rotating frame. This is an approxi- 


mation to the physical problem in which the injector 
flow upstream of the rotor is actually rotating with re- 
spective to the rotor in rotating frame. The correct mod- 
eling of the physical problem would require unsteady 
time-accurate simulations using a sliding interface to 
account for the interaction between the injector and ro- 
tor, or solve an isolated rotor with an unsteady inflow 
boundary condition. The purpose of this work is to 
verify the unstructured flow solver to model flow injec- 
tions with the source term method, and to predict (or es- 
timate) the impact of tip injections on the rotor perfor- 
mance, rather than to provide a detailed analysis of the 
unsteady interaction between the injector and the rotor. 
The computational mesh has about three million points 
for both injector and rotor, and is partitioned into 96 
blocks for parallel computations. 



Figure 9 Three Blade Passage Grid for Rotor 35 with 
Tip Injector 


Rotor 36 Oporating Chanctarwtic Corvw 



Figure 10 Predicted Rotor 35 Operating Ranges at 100% 
Design Speed 


Figure 10 shows the predicted rotor performance at 
the design speed with and without the tip injection. The 
total massflow shown in Figure 10 is the sum of the inlet 
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corrected massflow and the injector corrected mass- 
flow, which is calculated based on the following formu- 
la 


J T inlet/ T standard ( 
m total M inlet 77 > ' 

' inlet/ * standard 


JT injector /^ j 


m injector p t p 


standard 


( 8 ) 


where T is the total temperature and P is the total pres- 
sure. Computed choke massflows with and without tip 
injection seem to be smaller than the measured values. 
The reason for this may be due to the pressure loss cased 
by the extended shroud and hub surfaces upstream of the 
rotor. The total pressure and the total temperature need 
to be adjusted at the inlet to match the correct state for 
the measurement. Computed baseline stall massflow 
(no injection) is 19.18 kg/s. Computed stall massflow 
with the flow injection is extended to 17.38 kg/sec. The 
stall back pressure (on shroud) is increased from 
133,850 Pa in the case of no injection to 138,676 Pa in 
the case with injection. The stabilizing effect of the tip 
injection on rotor performance is clearly demonstrated. 
The predicted extension of operating ranges for Rotor 
35 is very similar to the numerical results by Suder, et 
al. [5] using the APNASA code, and also matches with 
the experimental trend. It is noticed that no total pres- 
sure loss occurred when the flow was injected through 
the casing in the axial compressor rotor. This is different 
from what was observed in the experiment for a high- 
speed centrifugal compressor, where a total pressure 
loss occurred when the flow was injected in the diffuser 
[ 20 ]. 



Figure 1 1 Static Pressure Contours at 98% Span of Rotor 



Figures 11 and 12 show predicted static pressure 
and axial velocity contours with tip injection at a loca- 
tion of 98 percent span of the rotor. The interaction be- 
tween the injected flow and the primary flow is clearly 
observed. The dash line indicates the interface between 
the injector grid and the rotor grid. It should be noticed 
that the flow field before the interface is viewed in the 
fixed frame, and the flow field after the interface is 
viewed in the rotating frame. The injected flow is 
pointed to the blade passage direction, rather than to the 
axial direction when it is viewed in the rotating frame. 

The effect of tip injections in suppressing the stall in- 
ception in the rotor is illustrated in Figure 13. The pre- 
dicted axial velocity contours with and without tip in- 
jections at the same back pressure (133,850 Pa on the 
shroud) are shown in Figure 13 (a) and (b). Low veloc- 
ity region occupied almost half of the blade passage at 
this back pressure in the case without tip injection, 
which led to continuous reduction of the massflow, and 


Figure 12 Axial Velocity Contours at 98% Span of Rotor 

eventually to the stall condition. The flow field was well 
retained as the flow was injected upstream of the rotor, 
which increased the momentum of flows in the tip clear- 
ance region. A stable flow condition was maintained 
even at a higher back pressure, and the stall point was 
effectively extended. 

As mentioned before, the current simulations only 
provide an estimate of the impact of tip injections on the 
performance of an axial compressor rotor. It is believed 
that the unsteady effect between the injector and rotor 
plays certain role in the stall inception process. Further 
more, because stall cells could rotate around the annulus 
of the compressor in the stall inception process which 
invalids the assumption of spacial and temporal period- 
icity, the correct way to model the stall inception pro- 
cess requires full annulus simulations of the compressor 
system. However, the current work demonstrated the 
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(a) Axial Velocity Contours Without Tip Injection 



Figure 13 Axial Velocity Contours on Cutting Plane 
Through Rotor Passage at a Back Pressure of 
133,850 Pa 

capability of the unstructured grid flow solver to simu- 
late the flow injection in a high-speed axial compressor 
rotor, which is useful for future investigations on stall 
inception process and stall control techniques for a full 
annulus simulations of compressor systems. 

Conclusions 

Numerical methods to model tip injection flows in 
an isolated compressor rotor were presented based on a 
Reynolds averaged Navier-Stokes arbitrary Mach num- 
ber flow solver. Simulations for a high-speed transonic 
compressor rotor, NASA Rotor 35 with and without tip 
injections were assessed with the experimental data. 
Computed operating ranges of Rotor 35 at 80 and 100 
percent of the design speed matched with the experi- 
mental trend. The stability enhancement of tip injec- 


tions on the rotor performance was demonstrated in the 
current work. Although the correct modeling of stall in- 
ception process requires time-accurate full annulus 
simulations of the compressor system, the methodology 
and numerical solutions presented in this work can be 
instructive for future investigations and optimizations 
of stall control techniques in gas turbine compressor 
systems. 
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Section II 


Sensitivity Analysis Enhancements and Results with the 
U 2 NCLE Software 



Sensitivity Derivatives 

For aerodynamic sensitivity analysis, the state equation is a system of nonlinear partial differential 
equations (PDE) expressing the conservation of mass, momentum, and energy. Differentiation of this 
system of PDE (i.e., sensitivity analysis) can be performed at one of two levels. In the first method, 
termed the continuous or variational approach, the PDE are differentiated prior to discretization, either 
directly or by introducing Lagrange multipliers which are defined as a set of continuous linear equations 
adjoint to the governing PDE. Subsequently, these directly differentiated or adjoint equations are 
discretized and solved. In the second method, termed the discrete approach, the PDE are differentiated 
after discretization. The discrete approach may also be cast in either a direct or an adjoint formulation, 
and the reader is referred to Hou et al. for a comprehensive presentation of both discrete formulations. 
For a more detailed discussion on the continuous approaches to aerodynamic design optimization the 
interested reader is directed, for example, to Refs. 2-4 
In general, the objective function and constraints may be expressed as Fj{p k ,X , q). Here, Q is the 

disciplinary state vector that is produced from the CFS, X is the computational mesh over which the 
PDE are discretized, and P* the vector of design variables. The sensitivity derivatives of these functions, 


VFj , may be expressed as 


dFj 

VF; = — L 

dp* 


'dFj'' 

dx 

b 

r dF j > 

U* J 

dP* 

UeJ 


dP* 


7 = 1, NCON + 1 


( 1 ) 


To compute the sensitivity derivatives in Eq.(l), the sensitivity of the state vector, d<2/dp* , and the grid 
sensitivity terms, dx/dp* , are needed. This approach to sensitivity analysis is referred to as the direct 
differentiation method. The number of linear systems needing to be solved is equal to the number of 
design or independent variables, NDV. If in the design problem under consideration, the sum of the 
objective function and constraints is less than the number of design variables (i.e., NCON+1 < NDV), a 
more efficient alternative approach may be formulated. This method is referred to as the adjoint variable 
approach, and may be written as 


VFj 



(dFj) 

a* -xL 

’ dR dR dX ' 
1 — 

dP* 

[ax J 

F ‘ 

dp* dx ap* 


7=1, NCON + 1 


( 2 ) 


where R represents the disciplinary state equation, and X F . are adjoint vectors defined in such a way as 


to eliminate the dependence of the objective function and constraints on the sensitivity of the state 
vector. This method requires the solution of NCON+1 linear systems. Hence, when the optimization 
problem formulation results in NDV > NCON+1, the adjoint variable approach should be used, whereas 
NDV < NCON+1 the direct approach. Methods used to perform the sensitivity analysis are discussed in 
subsequent sections. 


Discrete ( Semi-Analytic ) Approaches 

For the discrete-direct approach (Eq.l), the sensitivity of the field variables are required, and for the 
discrete-adjoint approach (Eq. 2) the adjoint vectors are needed. To obtain these, the discrete residual 
vector for a steady-state solution may be recast as 

F(p*,x(p*),£>(P*))=o (3) 


where the explicit and implicit dependencies of the residual on the state vector Q(fi k ), the computational 
mesh X (P* ) , and the design variables P* are asserted. 

In the discrete-direct approach, Eq.(3) is directly differentiated with respect to the design variables to 
produce the following linear equation 


or, rearranging 


dR_ = dR_ dR_dX_ dR dQ 
d$k dP* dX dP* dQ dp* 

dR dQ _ f dR dR dX ' 
dQ dP* ~ [ap* + dX dp* 


(4a) 

(4b) 


where d R/dQ and dR/dX are the Jacobian matrices evaluated with a converged (steady-state) solution. 
The solution of Eq.(4b) poses the difficulty of solving a large linear system of equations for each design 
variable. Furthermore, due to the fact that this is a linear system, the linearizations of the residual with 
respect to the state vector and the grid (i.e., the Jacobian matrices) must be exact. The detriments of 
using inexact linearizations in Eq.(4b) have been explored in Ref. 5. Solving these systems, however, is 
made more tractable when the above equations are recast into what has been termed the incremental 
iterative form [6,7] as follows 


(dQ) 


dR dR dX dR dQ 

V H” 

^ap* j 
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where A may be any convenient approximation to the higher-order Jacobian which converges the linear 
system. This is because the equations are now cast in delta form, with the physics contained in the right- 
hand-side vector. It has been found that the first-order Jacobian works well for use in the coefficient 
matrix of Eq.(5a); most CFS codes also use the first-order Jacobian for this purpose. 

Two particularly attractive features of the incremental iterative strategy are that (i) a more diagonally 
dominant matrix may be used to drive the solution of the linear systems (as opposed to the sometimes 
ill-conditioned higher-order Jacobian), and (ii) the higher-order Jacobian now resides on the right-hand- 
side of the equations and may be dealt with in an explicit manner. When in this form, only the fc-vectors 
resulting from the matrix-vector product of [dR/dQ:Jf)Q/d ^ k ) are of concern. Hence, CPU time and 
memory efficient methods for constructing the exact matrix-vector product can be utilized. To this end, 
higher-order spatially accurate discrete-direct sensitivity analysis procedures have been developed and 
reported in Ref. 8. 

The discrete-adjoint variable formulation begins by combining Eq.(4a) from the direct differentiation 
method with the sensitivity derivatives in Eq.(l). From this, adjoint vectors may be conveniently defined 
such that the sensitivity of the field variables are no longer needed. Omitting details, the end result 
requires the solution of the following linear systems for the adjoint vectors 

dFj Y 

— f j = 1, NCON + 1 (6) 

x < d QJ 

The adjoint variable approach may also be recast in incremental iterative form, however, for brevity 
these details shall be omitted. It should be noted that all the same linearizations required by the direct 
approach are required by the adjoint variable method; they are simply transposed and used at different 



stages in the computation of the sensitivity derivatives. Hence, a sensitivity analysis code which was 
developed for either of the discrete approaches may be modified to produce the other. 

Numerical Approaches 

The task of constructing exactly or analytically all of the required linearizations and derivatives by 
hand for either the discrete-direct or discrete-adjoint approach and then building the software for 
evaluating these terms can be extremely tedious. This problem is compounded by the inclusion of even 
the most elementary turbulence model (for viscous flow) or the use of a sophisticated grid generation 
package for adapting (or regenerating) the computational mesh to the latest design. One solution to this 
problem has been found in the use of a technique known as automatic differentiation. Application of this 
technique to an existing source code, that evaluates output functions, automatically generates another 
source code that evaluates both output functions and derivatives of those functions with respect to 
specified code input or internal parameters. A pre-compiler software tool, called ADEFOR (Automatic 
Differentiation of FORtran, Bischof et al. 9 ), has been developed and utilized with much success to 
obtain complicated derivatives from advanced CFD and grid generation codes [10,11]. The use of 
ADIFOR produces code that, when executed, evaluates these derivatives of the output functions via a 
discrete-direct approach, referred to as forward-mode automatic differentiation. More recently, 
automatic differentiation software has emerged that enables the derivatives to be evaluated with a 
discrete-adjoint approach [12,13]. This type of automatic differentiation is known as reverse-mode. A 
detailed and concise review on the use of sensitivity analysis in aerodynamic shape optimization has 
been reported by Newman et al. 14 , the reader is directed to this source for a discussion on the methods 
presented above. 

The methods discussed above require differentiation of the CFS software, either by hand or with pre- 
compiler software. Other methods to obtain sensitivity derivatives are based on numerical techniques. 
The simplest numerical technique is the finite-difference approximation. Another is a new numerical 
technique, pioneered at MSU [15] for performing disciplinary and multidisciplinary sensitivity analysis 
which uses complex variables to approximate derivatives of real functions. This method is based on 
ideas that where explored over three decades ago by Lyness and Moler 16 and Lyness 17 , and recently 
revisited by Squire and Trapp 18 . Both numerical approaches will be discussed below. 

For a central finite-difference approximation to the derivative, one may expand the function in a 
Taylor series about a given point using a forward and a backward step, and then subtracting to yield 

dF [F(x + ft)-.F(.x-ft)] ft 2 d^F h^d^F 

dx~ 2 h 3/ dx 2 5/ dx 5 { 

This expression for the derivative has a truncation error of o(/z 2 ). The advantage of the finite-difference 
approximation to obtain sensitivity derivatives is that any existing code may be used without 
modification. The disadvantages of this method are the computational time required and the possible 
inaccuracy of the derivatives. The former is due to the fact that for every derivative, two well-converged 
CFS solutions for the function evaluations are required (i.e., 2 *NDV analyses). In the case of nonlinear 
fluid flow, for example, these solutions may become extremely expensive. The latter disadvantage is 
attributed to the sensitivity of the derivatives brought about by the choice of step size. To minimize the 
truncation error one selects a smaller step size, however, an exceedingly small step size may produce 
significant subtractive cancellation errors. The optimal choice for the step size is not known a priori, and 
may vary from one function to another, and from one design variable to the next. 

Instead, if the function is expanded in a Taylor series using a complex step as 



(8a) 
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Solving this expression for the imaginary part of the function yields 

dF _ Im[F(x + /»)] ft 2 d 3 F ft 4 d 5 F 
h + 3/ ^3 5/ ^5 + 


(8b) 


This expression for the derivative also has a truncation error of d(/z 2 ). By evaluating the function with a 
complex argument, both the function and its derivative are obtained, without subtractive terms, and thus 
cancellation errors are avoided. The real part is the function value to second order. This technique has 
been labeled as the complex Taylor series expansion (CTSE) method [15]. 

Previously, a disadvantage of the complex variable approximation was the increased runtime required 
by the evaluating routines when run with complex arguments. With current compiler options, and when 
the entire CFS code was converted to complex, that run time was on the order of three times the cost of 
the original solver. Automatic differentiated versions of analysis software, using say ADIFOR, incur 
about the same time penalties. Currently, the run time using the complex variable approximation costs 
little more than hand differentiated, direct-discrete approaches. A detailed comparison of sensitivity 
derivatives from automatic differentiated software and using the complex variable technique can be 
found in Ref. 19. The reader is directed to this source for a more complete discussion on the competing 
methods. 

The advantages of the complex variable approximation for obtaining the derivatives are numerous. 
First, like the finite-difference approximation to the derivatives, very little modification to the software 
is required. All the original features and capabilities of the software are retained. Thus, user experience 
with the current software is not lost, and ongoing advancements and enhancements can be readily 
introduced into subsequent versions without extensive modifications or re-differentiation. This is in 
direct contrast to hand or automatically differentiated codes where any modification to the original 
software will require re-differentiation. This advantage is extremely useful in the problem formulation 
stages of the design process when new objective functions and constraints are being explored. Second, 
this method has been shown in Ref. 20 to be equivalent to a discrete-direct approach, either from 
automatic differentiation or hand differentiated codes solved in incremental iterative form, in the way 
that the state vector and its derivatives are being solved for simultaneously. When solving the state 
equation, the state vector resides in the real part and the derivatives in the imaginary part. Hence, fully 
converged flow solutions are not required to obtain derivatives of sufficient accuracy for design. Finally, 
the complex variable approximation to the derivative is not sensitive to the step size selection and only 
requires step sizes that avoid excessive truncation error; thus, this method demonstrates true second- 
order accuracy [15,21]. The CTSE method may also be used to construct any linearization or Jacobian 
required for CFS [22]. This may be particularly useful when complicated flux functions or turbulence 
models are being developed. In addition, the complex variable technique can be used to compute second 
derivative information using available data [19]; however, these computations are subject to cancellation 
errors. v - 


Sensitivity Derivative Validation Study 

In this section, validation of sensitivity derivative calculations for a three-dimensional wing are 
presented. Sensitivity derivatives with respect to lift and drag coefficients are computed with both the 
discrete-direct and discrete-adjoint formulations for low- and high-Reynolds number viscous, laminar 
flows. The configuration chosen is the standard benchmark ONERA M6 [23] wing, and the flow 


conditions assumed are a Mach number of 0.52 with a 2.0° angle-of-attack. The y-coordinate of a grid 
point on the wing surface is selected as the independent variable. 

In order to validate the computed sensitivity derivatives with finite-difference values, a detailed step- 
size study must first be performed using finite-differences. This is required in order to find an acceptable 
step-size for finite-differences, which are subjected to both truncation and cancellation errors. This study 
for sensitivity of lift and drag coefficients, at Reynolds numbers of 100 and 510 6 , for forward, 
backward, and central finite differences are shown in Fig. 1. Step-sizes above and below that shown are 
subjected to the aforementioned errors, and can not be presented on the scale chosen in these plots. 
Selecting a finite-difference step-size of 10' 8 to compare with, the computed sensitivity of lift and drag 
for the low- and high-Reynolds number simulations are given in Table 1 and Table 2, respectively. 
Observe that the accuracy of the computed derivatives, using both the direct and adjoint formulations, 
agree to within 1% of the central finite-difference value. Although this accuracy is exceptional, it is 
more than that which is required to perform computational design, error estimation, or uncertainty 
analysis. Note, however, that obtaining accurate sensitivity derivatives using finite-differences is 
usually difficult, if not impossible. Hence, standard test configurations such as the ONERA M6, which 
have been used for code validation and numerous sensitivity analysis studies, are a prerequisite. 

The direct approach to sensitivity analysis produces the derivative of state vector with respect to the 
independent variable. Hence, it physically represents the rate of change of dependent to independent 
variable. The adjoint vector, on the other had, does not have an obvious physical interpretation, but it 
does represent the rate of change in the objective/constraint function (see Eq. 6) to the state-equation. 
Thus, the adjoint vector has utility in error estimation, error correction, and grid adaptation. The use of 
the adjoint vector for this purpose has been explored in Ref. 24. 

Conclusion 

In the current research, both the discrete-direct and the discrete-adjoint formulations for sensitivity 
analysis have been incorporated into the U 2 NCLE software and validated for compressible viscous, 
laminar flows on a benchmark configuration. This validation study included a detailed finite-difference 
step-size study in order to compare with accurate numerical derivatives. It should be noted that for 
rotating machinery the rotational speed, and therefore the dynamic grid terms, are not a function of the 
geometric shape of the configuration. Hence, the linearizations for the discrete sensitivity analysis 
performed in this study are valid for that type of analysis as well. 
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(b) Sensitivity of drag, Re=100. 
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(d) Sensitivity of drag, Re=510 6 . 


Figure 1: Step-size study for finite-difference sensitivity derivatives. 


Table 1: Comparison of sensitivity derivatives, Reynolds number of 100. 


Sensitivity 

Direct Analysis 

Adjoint Analysis 

Finite Difference* 

dC L / 
/dP 

0.000565 

0.000565 

0.000565 

dC D / 

/dfi 

0.000927 

0.000926 

N -v 

0.000922 


* Central finite difference for step size of 10' 8 



Table 2: Comparison of sensitivity derivatives, Reynolds number of 5T0 6 . 


Sensitivity 

Direct Analysis 

Adjoint Analysis 

Finite Difference* 

dC L / 

/dfi 

0.0000691 

0.0000697 

0.0000693 

dC °/„ 

/dp 

-0.0000102 

-0.0000102 

-0.0000102 


* Central finite difference for step size of 10' 8 
















